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We present numerically exact energy estimates for two-dimensional electrons in a parabolic con- 
finement. By application of an extension of the recently introduced many-body diffusion algorithm, 
the ground-state energies are simulated very efficiently. The new algorithm relies on partial an- 
tisymmetrization under permutation of particle coordinates. A comparison is made with earlier 
theoretical results for that system. 
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I. INTRODUCTION 



In this paper we apply the recently developed quantum Monte Carlo algorithms to an important model 
for interacting charge carriers in quantum dots. The fabrication of novel miniaturized semiconductor structures on 
nanometer scale has shed light on a variety of advanced physical systems and devices in which the classical description 
' of electronic properties breaks down. The band structure of quantum wells, multiple quantum wells and superlattices 
H makes mobile carriers locate parallel to the semiconductor interfaces and hence induces a quasi-two-dimensional 
confinement H . Electrons in sandwiched semiconductor layer structures can be confined perpendicular to the growth 
direction If the lateral confinement is of the same order as the electron wavelength, the electrons have essentially 

no free direction left. The resulting quasi-zero-dimensional structure is addressed as a quantum dot or an artificial 
atom H . Mostly, the confinement effect along the growth direction of the layers is much stronger than perpendicular 
to it. Accordingly, on the microscopic scale, the electrons form two-dimensional disk-shaped objects 0. Step-like 
potential structures, being induced by steps in the conduction band edge, became a popular object of study in 
the framework of self-assembled quantum dots ||. In field-effect confined quantum dots, experimental far- infrared 
transmission spectra indicate the occurrence of parabolic-shaped confinement potentials ||. Depending on the gate 
voltage applied, the number of electrons may discretely vary from zero to a large number of electrons ]lO|JTlJ . 
O ■ Motivated by the fabrication techniques of quantum dots, we here focus on a commonly used model of N interacting 
\ electrons in a harmonic confinement potential, 

•i-H . 



H=^ + V c (f) + V int (f) (1) 
where 

N-l N 2 
i— 1 j—i J 

Here, u denotes the frequency of the harmonic confinement potential of the dot, fi the effective mass of the electrons 
and e the dielectric constant of the material. For convenience, we set h = u) = e = p,= l. The only parameter which 
then enters in (1) is 7 = e _1 . It serves as a parameter to adjust the strength of the electron-electron repulsion for a 
given dot radius. 

The derivation of the energy eigenfunctions and eigenvalues of (|l|) turns out to be challenging in cases where many 
particles are involved. Even in the classical case, solutions remain elusive, and the use of analytical approximations 
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or computational methods is indispensable to reliably predict the energy spectrum of the system. The consideration 
of identical particles makes analytical or numerical treatments still more difficult. Much effort with different methods 
has been pursued to investigate the physical properties of quantum dots 12 ljj. Various approaches have been 
performed using the Thomas- Fermi approximation |fj"4| , many-body perturbation theory |fj"5fl , or the Hartree- and 
Hartree-Fock Ansatz |f^] . Bearing in mind that exact analytical approaches are generally restricted either to the limit 
of non-interacting quantum- mechanical particles or exclusive cases |17|| , some methods, as e.g. the Pade approximants 
jl8],[L9| or the renormalization perturbation theory [^0 21 1, attempt to connect these limits. Exact numerical solutions 
are possible by diagonalization of the Coulomb interaction 15 Tqj23-p3. Due to their challenging numerical expense, 



diagonalization techniques are limited to few-electron systems. In particular for many-particle systems, for which 
diagonalization techniques quickly exceed the capacities of today's computers, Monte Carlo techniques turn out to be 
very useful ]25| ]. Monte Carlo approaches have been applied to electrons in a parabolic potentials in refs. fp6|-p8[. Apart 
from their simplicity and flexibility, the power of Monte Carlo methods lies in the possibility to obtain estimates which 
converge to the exact values with a known statistical error. Although the computational effort needed for Monte Carlo 
techniques to obtain a given accuracy does not necessarily grow exponentially with the particle number, it mostly does 
if identical fermions are involved. As a consequence, exact fermion Monte Carlo methods are limited to few-particle 
systems. In what follows, we apply a Monte Carlo method which scales favorably and is nevertheless exact. 

The paper is organized as follows. In section II we will discuss the principle of many-body diffusion (MBDF) 
po| p2|. This formalism is applied to derive random processes to numerically predict some energy eigenvalues of the 
model Hamiltonian. A feasible implementation scheme is discussed. Section III gives the outcome of our numerical 
analysis, and a comparison to other work is made. Finally, section IV concludes the present article. 



II. THE MANY-BODY DIFFUSION ALGORITHM 



The model is approached with a recently developed quantum Monte Carlo method, the many-body diffusion al- 
gorithm (MBDA) §],§. The MBDA is based on the many-body diffusion formalism (MBDF) jH^. Here we 
will restrict ourselves to a brief outline of the underlying concepts. In the MBDF, the propagator of N interacting 
identical particles in d spatial dimensions is written as a Feynman-Kac functional |33| over a symmetrized process, 
i.e., as a Euclidean-time path integral over the diffusion process of N identical free particles with superimposed 
potential-dependent exponential weights. For a given irreducible symmetry representation S, the propagator 



K s (rf,r;ri) = E- 



(3) 



is hence represented as an average over all paths starting in as indicated by the averaging index E- t , and ending a 
Euclidean time lapse t later in Tf, as denoted by the indicator I(ji^ = - f y The symmetry representation S determines 

the construction principle for the underlying <i/V-dimensional diffusion process ^R S (r) ; r > o\. 

In the MBDF, a detailed analysis has been performed of the diffusion process of free identical particles, and 
the role of the potential symmetry has been pointed out. It was found that for coordinate-symmetric potentials, i.e., 
potentials invariant under the permutation of the Cartesian particle coordinates, and for certain irreducible symmetry 
representations S, the total propagator separates into a sum of stochastically independent sub- propagators. The 
importance of this coordinate-symmetry shows two- fold. First, it allows to easily generalize the diffusion process of 
free identical particles to the process of interacting identical particles. Second, dealing with identical fermions, the 
sign problem is strictly avoided by numerical procedures based on exclusively positive walkers. The Feynman-Kac 
formulation (^) indicates the relevance of the free diffusion process. For any coordinate-symmetric potential V(f), 
the symmetry properties of the total diffusion process correspond to those of the free one. Rather than interfering 
with the role of the potential, different irreducible symmetry representations S do completely determine the structure 
of the free diffusion process. Correspondingly, one may in principle approach various eigenstates of the system by 
formulating the free diffusion process in the appropriate symmetry representations. 



In 1 31 1, the free density matrix of N identical particles is decomposed into corresponding one-dimensional ./V-particle 
density matrices. In one dimension (ID), the introduction of an ordered iV-particle state space = {x\ > > 
■ ■■ > xjv} projects the density matrix on a mathematically well-defined expression for both bosons and fermions. 
On that basis, the MBDF introduces the fermion diffusion process Xf and the boson diffusion process Xf, as a 
Brownian motion on the irreducible state space Dn with absorbing respectively reflecting boundary conditions. The 
processes Xf and Xb serve as key ingredients for the multi-dimensional formulation. With the decomposition into 
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one-dimensional fermion and boson diffusion processes, a scheme has been introduced to sample the free density 
matrix for specific symmetric and antisymmetric symmetry representations sign-problcm-free. 

The symmetry constraints specified in [[l 31 rely on the complete (anti-)symmetrization along the Cartesian co- 
ordinates. This scheme addresses a particular excited fermion state. To illustrate this idea, consider N identical 
free fermions or bosons in two dimensions (2D) with unit mass, for which the density matrix can be expressed as a 
determinant or a permanent: 



Pt(rf,r;ri) 



det 

j,k=l,N 



p(dr/,r;r i fe ) resp. p h (77, r; n) 



perm 

j,k=l,N 



Pd 



( r 7 



■r k 



where 



Pd 



2ttt 



exp 



(f/ 



■k\2 



2t 



Complete particle (anti)symmetrization along the Cartesian directions then leads to the following representation for 
N two-dimensional identical fermions 



Pi(rf,T;ri) = det 

j,k=l,N 



Pd 



pt (x f ,T-,Xi)p h (vf,T;yi) + p h (x f ,T;xi)p { (jZ/.r;^) 



Applying an analogous decomposition for three spatial dimensions, excited state energies of up to 20 harmonically 
interacting identical fermion oscillators have been efficiently simulated within a statistical accuracy of about 0.1 
percent W . Quadratic particle interaction has also been considered in the framework of quantum dots fl35| but is not 
studied here. By an extension of the many-body diffusion principle, it is also possible to extract the corresponding 
ground-state energies. Again, particle (anti)symmetrization along the Cartesian directions plays the central role in 
that formulation. However, in contrast to fljjn}, not all particles are necessarily (anti)symmetrized in the same 
direction. A detailed discussion of the underlying formalism is beyond the scope of the present paper and will be 
published elsewhere. 

The derivation of the sampled functional separates into two parts, the decomposition of the corresponding free 
density matrix and symmetry considerations on the potential. The underlying free diffusion principle is best explained 
for the example of three identical free fermions in two dimensions |34j] . In this case, the infinite- time limit of the free 
density matrix reads 



lim 

r — >oo 



Pf 0"/,' 



Pd (r/, 



\e r f 



1 xi yi 1 x\ 

1 X 2 V2 1 X' 2 

1 x 3 y 3 1 x' 3 

( [{xi -x 2 )(y 2 - 
f [(x 2 - xa)(yi 
- [(xi - x 2 )(y 2 
\ - [(x 2 - x 3 )(yi 



y'2 



0(r~ 3 ) 

2/3 ' 

ife)] [(x[-x' 2 M-y' 3 )] 
-V2)] [(x' 2 -x' 3 )(y[-y' 2 )] 
-2/3)] [(x' 2 -x' 3 )(y[-y 2 )} 

- V2 )] [K-4)(^-y 3 )] 



+ 0(r- 3 ). 



(4) 



As our concern is the generation of equilibrated samples to derive properties of the (lowest available) eigenstate, the 
use of this asymptotic limit is justified. Emphasis in our approach -as mostly in random- walk or diffusion Monte Carlo 
approaches- is thus on the long-term distribution rather than on the equilibration process itself. Eq. (^) involves 
two interdependent stochastic processes. They both include a Brownian motion, but distinguished by their respective 
positive domains D\ and D 2 , 

Di = {xi > x 2 ; y 2 > y 3 } and D 2 = {x 2 > x 3 ; y x > y 2 } , 



and boundary conditions. Apart from an adapted distinguishable-particle diffusion, a jump process must be realized 
to take into account the process interdependencies g]. In practice, during an evolution cycle, a walker associated 
to a particular type of process might be assigned to the other type of process according to the locally dependent 
process transition rates. In the present case, in which we are interested in the derivation of the ground-state energy, 
single-process evolution is sufficient for our needs. Indeed, as long as the time decay rates of the individual processes 
are identical, the energy eigenvalue predicted by any of the single processes is the same as the one predicted from 
their combination. 

Table 1 shows the symmetry configurations for the numerical simulation of the ground-state energy of three and 
six two-dimensional non- interacting spin-polarized harmonic fermions. The ground-state energy of three identical 
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fermions in 2D, e.g., can be simulated by anti-symmetrization of pairs of Cartesian coordinates, namely {x\,X2) and 
(2/2,2/3)- Analogously, one might conglomerate the coordinates (2:2, £3) and (2/1,2/2)- For six identical fermions in 2D, 
e.g., one simultaneously anti-symmetrizes (x±, X2, £3), (xi,Xe), (2/3,2/4,2/5) , and (2/2,2/6), and so forth. 

The advantage of the (generalized) many-body diffusion approach lies in the efficient sampling of the required 
probability densities. As mentioned above, symmetry and anti-symmetry in one spatial dimension can be achieved by 
the definition of reflecting and absorbing boundary conditions on a distinguishable-particle Brownian motion. Neither 
do we have to deal with negative transition amplitudes, nor is our approach slowed down by sampling determinants. 
A detailed description of the algorithmic realization has been reported in jl) . The Euclidean-time evolution according 
to Kg (r /,T;f,) is simulated in sufficiently small time steps e by the repeated application of the following two- 
step procedure: a) given (2A r -dimensional) initial system configurations sample final ones f>(r + e) = fi(r) + 
Sr(e), where Sr(e) is randomly drawn according to the free identical-particle propagator, and b) apply the potential- 
dependent weights exp 

the weights exp 



-C<v(R S (<))d< 



randomly in a branching and killing procedure. The determination of 

— JJ +€ V (^R S (s)J in principle requires infinitesimal time steps e — > 0. Due to limited computer 

performance, however, this procedure is not practical, and reliable approximation schemes must be provided. In 
the present case of harmonic confinement and repulsive Coulomb interaction, the use of the Suzuki- Trotter weights 



exp^-i 



V 



R S ( T )^J + V (^R S (r + e)^ I is satisfactory for realistic time steps of O.OOf /Hartree (H). The essential 
requirement for the efficient approach of many-body systems with the outlined free diffusion construction principles is 
the coordinate-symmetry of the potential involved. This condition holds for both the confinement and the interaction 
potential of our model 

V(i, j e {1,2,3} V {4,5,6}) : V int (f) = V int (P^f) = V int (P^f) A V c (f) = V c (P^r) = V c (P^f). (5) 

In (||), the operators P^ and P*' J interchange the ith and the jth x- respectively y-coordinates. 

It should be emphasized that this invariance of the Hamiltonian under the interchange of the x and y coordinates of 
any two particles no longer applies if the confinement potential is replaced for instance by a Coulomb potential. The 
potential then generates transition rates between different types of walkers, as discussed in detail in Although the 
principle of sign-problem free diffusion remains valid for such systems, the detailed analytical analysis of the process 
interdependencies becomes cumbersome. With our present approach of Carthesian decomposition, more than six 
electrons become almost intractable in practice. New purely numerical techniques are currently under development to 
perform the required symmetry decompositions, avoiding the tedious and unpractical analytical bookkeeping of the 
Carthesian decompositions for the different types of walkers. Preliminary studies reveal that this numerical analysis 
requires in general the evaluation of N x N determinants relating the initial and final positions of the walkers. This 
N cost in computation time is presumably overcompensated by a factor 1/N\ due to the reduction of the state space. 
Within the present status of our approach however, closed shell systems with potentials satisfying the symmetry 
condition (||) are tractable for as many as 20 electrons. 

III. RESULTS AND DISCUSSION 

Table 2 gives the ground-state energy estimates predicted for closed-shell systems of six, twelve and twenty unpolar- 
ized electrons as a function of 7 1 / 3 . The energies obtained are supposed to be exact within the numerically estimated 
statistical error. In the limit of zero electron Coulomb repulsion, 7 — > 0, (1) is identical to a system of non-interacting 
fermion oscillators. The arrangement of single-particle harmonic oscillator solutions into a Slater determinant then 
induces the energy limits E^ {) = 10 H, Eq^° = 28 H and Eq^° — 60 H for six, twelve and twenty electrons, respec- 
tively, whereas the neglect of spin statistics would yield Ej^ = N H, with N indicating the particle number. In the 
opposite limit, for infinitely large electron-electron repulsion, 7 — > 00, the electrons behave classically and one expects 



them to arrange in the form of a Wigner lattice |36|. With increasing 7, the average electron-electron distance grows, 
and the influence of spin-statistics weakens. Accordingly, as 7 grows, the energy gap between different excited states 
is expected to decrease substantially. This physical behavior is recovered by our numerical data. The limit of zero 
electron-electron repulsion is accurately simulated, and a smooth transition to high 7 is found. A comparison of the 
numerical energy eigenvalues and the energy of distinguishable particles (see Table 3) indeed indicates the irrelevance 
of quantum statistics as 7 — > 00. 

It proves instructive to compare our energy estimates for the ground state of the unpolarized electron systems with 
the Pade approximants reported in ref. p9[ . The relative deviation of the Pade approximants with respect to our 
numerical estimates Eq for the case of six unpolarized electrons are also studied. For both zero and very large 7, 
the Pade approximates ^-3,2(7) and ^4,3(7) match our prediction. For intermediate 7, the Pade approximants ^4,3(7) 
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for the ground state energy introduce a systematic error of up to almost 4 per cent. A comparison of our numerical 
estimates for the twelve and the twenty-particle system with the corresponding Pade approximates [fl9f indicates 
an analogous qualitative picture. The extreme relative deviations are of the same order. The lack of an accurate 
description of the system (1) for intermediate regions of the electron-electron- repulsion parameter 7 is typical for a 
variety of analytical approximations. 



IV. CONCLUSIONS 



In the present paper, we apply a generalization of the recently reported many-body diffusion formalism to a system 
of interacting electrons in a parabolic confinement. The method is illustrated explicitly for closed-shell configurations 
of six, twelve and twenty unpolarized electrons for which the ground-state energy is numerically predicted. The 
algorithm proceeds without the use of analytical approximations. Apart from a small but controllable systematic 
error, the energy values are numerically exact within a computed statistical error of a few per mil. The feasibility 
of our approach is indicated by the comparison of our ground-state energy estimates with the corresponding Pade 
approximants calculated by A. Gonzales et al. JljJ. Analogous results |37) were obtained by another Monte Carlo 



technique [E8 38 that strongly reduces the noise due to the sign problem. 



Regarding the methodological aspect, this work introduces an algorithm which allows to efficiently simulate the 
ground-state energy of closed-shell configurations of electrons exposed to coordinate-symmetric potentials. The scheme 
strictly avoids the fermion sign problem p^ffOt by the definition of a Brownian motion on a state space with the 
appropriate boundary conditions. The resulting random process can be realized by stable diffusion of purely positive 
walkers. 

The formulation of a sign-problem-free algorithm for the quantum-dot model (1) is remarkable, since almost all 
quantum Monte Carlo algorithms face massive difficulties with the general description of significantly correlated 
continuous quantum systems. The reason for this is the fermion sign problem, which generally thwarts reliable 
stochastic many-fermion treatments. Although potentially exact, the transient estimation of eigen-energies experiences 
a serious inefficiency due to an exponential decreasing signal-to-noise ratio. On the other hand, even the use of very 
accurate and efficiently implemented trial wave functions in diffusion or Green's function Monte Carlo variants |4l]] 
introduces considerable systematic errors due to the assumption of generally incorrect nodal surfaces. The MBDA 
avoids these problems by the definition of a diffusion process on a state space with absorbing and/or reflecting 
boundaries. Its derivation is based on a symmetry analysis of the Hamiltonian only. 
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G 



N 


antisymmetric constellation 


3 


(x 1 ,x 2 ); (2/2,1/3) 


6 


(x 1 ,x 2 ,x [i ); (x 4 x e> ); (1/3,2/4,2/5); (2/2,2/6 ) 


10 


(2/1,2/2,2/3,2/4); (2/5,2/8,2/9); (2/6,2/10) 

(£4, £5, £6, ^7); (2:3,^8,2:10); (£2,2:9) 



TABLE I. List of configurations used for the ground states of N identical 2D fermions. 



y/3 


0.2 


0.4 


0.6 


0.8 


1.0 


1.2 


1.4 


1.6 


£o(3T,3|) 
^o(6T,6|) 

£o(iot, 101) 


10.086(4) 
28.28(4) 
60.99(7) 


10.804(3) 
31.00(3) 
67.72(9) 


12.632(4) 
37.91(2) 
85.24(3) 


15.806(4) 
49.62(2) 
114.98(8) 


20.367(3) 
66.45(1) 
157.97(4) 


26.158(2) 
87.88(2) 
211.70(8) 


33.084(4) 
113.52(2) 
276.72(7) 


41.115(5) 
143.45(2) 
351.72(6) 



TABLE II. List of energy eigenvalues (a.u.) for model (1) obtained by the MBDA. 





£o(3T,3|) 


Sdist(3T, 3J.) 


£<>(6T,6j) 


S di st(6T,6|) 


£ (10T,10|) 


£ di st(10T, 104) 


7 1/3 = 2 


60.385(5) 


60.12(1) 


215.21(1) 


214.30(3) 


533.56(3) 


531.44(4) 



TABLE III. List of energy eigenvalues (a.u.) for model (1) obtained by the MBDA. 
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